Tunneling and delocalization effects in hydrogen bonded systems: a study in position 

and momentum space 



o 
o 

(N 

Ph. 
< 

(N 
o 

S 

> 

a 

'Til 

o 
o 



Joseph A. Morrone 
Department of Chemistry 
Princeton University 
Princeton, NJ 085iU 

Lin Lin 

Program in Applied and Computational Mathematics, 
Princeton Universtty 
Princeton, NJ 08544 

Roberto CaiQ 
Department of Chemistry and Department of Physics, 
Princeton University 
Princeton, NJ 08544 
(Dated: April 30, 2009) 

Novel experimental and computational studies have uncovered the proton momentum distribu- 
tion in hydrogen bonded systems. In this work, we utilize recently developed open path integral 
Car-Parrinello molecular dynamics methodology in order to study the momentum distribution in 
phases of high pressure ice. Some of these phases exhibit symmetric hydrogen bonds and quantum 
tunneling. We find that the symmetric hydrogen bonded phase possesses a narrowed momentum 
distribution as compared with a covalently bonded phase, in agreement with recent experimental 
findings. The signatures of tunneling that we observe are a narrowed distribution in the low-to- 
intermediate momentum region, with a tail that extends to match the result of the covalently bonded 
state. The transition to tunneling behavior shows similarity to features observed in recent experi- 
ments performed on confined water. We corroborate our ice simulations with a study of a particle 
in a model one-dimensional double well potential that mimics some of the effects observed in bulk 
simulations. The temperature dependence of the momentum distribution in the one-dimensional 
model allows for the differentiation between ground state and mixed state tunneling effects. 

PACS numbers: 71.15.Pd 
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I. INTRODUCTION 

The nature of the hydrogen bond plays a critical role in 
determining the behavior of biological and chemical sys- 
tems. The state of a proton that participates in hydrogen 
bonding, therefore, is a subject of great interest. This 
is often characterized in terms of position distributions 
such as the radial distribution function. The radial dis- 
tribution function quantifies the probability of one atom 
being located a certain distance from another atom in po- 
sition space. This quantity may be extracted from elastic 
neutron and x-ray scattering experiments (see e.g. Ref- 
erences [1133) ; s-nd may also be computed via molecular 
simulation^. 

Recently, the position space picture that is provided 
by the radial distribution function has been comple- 
mented by measurements of the proton in momentum 
space. Deep inelastic neutron Compton scattering ex- 



perimenti 



,5,6,7 



have uncovered this property in a variety 



of hydrogen bonded systems, including several phases of 
bulk water— water confined in nanomaterial o^^'^^ and 
biological systemsi^, ferroelectric o^*^'^^ and superprotonic 
conductorsi^. In these experiments, the observed mo- 
mentum distribution sharply deviates from classical be- 
havior, underlining the importance of nuclear quantum 



effects. These phenomena may be understood in terms of 
the intertwining of the momentum distribution with the 
potential energy surface that originates in quantum me- 
chanics due to the uncertanty relation between position 
and momentum. As an example of this property, consider 
the weakening of the oxygen-hydrogen covalent bond as 
the strength and stability of hydrogen bonding increases. 
Such effects are typically associated with the red-shift of 
the OH stretch in the infrared spectra, and may be ob- 
served in the shortening of the tail of the proton momen- 
tum distribution^. In certain systems, the proton may be 
shared equally between the donor and recipient oxygen, 
thereby forming a so-called "symmetric" hydrogen bond. 
If one follows the logic of the weakening covalent bond 
as the hydrogen bond strengthens, then one would ex- 
pect an increased narrowing of the tail of the momentum 
distribution as two weak "symmetric" bonds arc formed 
between the hydrogen and its neighboring oxygens. This 
is precisely the type of behavior that has been observed 
in experiments on water confined in carbon nanotubes^i 
and the superprotonic conductor Rb3lI(S04)2^. 

As nuclear quantum effects are essential in determin- 
ing the form of the momentum distribution, it would 
seem natural that tunneling phenomena may be stud- 
ied by means of this property. Indeed, neutron Comp- 
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ton scattering experiments have presented momentum 
distributions that report signatures of tunnehng in sys- 
tems such as the ferrolectric potassium dihydrogen phos- 
phate (KDP)^-, supercooled water—, water confined in 
sihca nanopores^^ and the hydration shell of globular 
proteinsi^. In these systems, a secondary maximum or 
shoulder is observed that has been related to a node in 
the tunneling direction of the momentum distribution. 
This feature is further supported by a simple analytical 
model of the ground state of the tunneling wavefunction^ . 
Furthermore, in the experiments that involve supercooled 
and confined wate r^i^^i^^ , an excess of kinetic energy is 
observed when compared to liquid water at ambient con- 
ditions. This is manifested in momentum distributions 
that possess higher density at larger momenta. 

The proton momentum distribution may also be ac- 
cessed via computer simulation. As discussed above, nu- 
clear quantum effects are essential to the computation of 
this property. Typical atomic simulations treat the nu- 
clei as classical point particles. In such simulations the 
momentum distribution is given by the Boltzmann dis- 
tribution and exhibits no dependence upon the potential 
energy surface. Nuclear quantum effects may be included 
within the Feynman path integral formulation of quan- 
tum statistical mechanics^^. In the discretized implemen- 
tation of this method, one quantum mechanical particle 
is mapped onto a number of classical replicas that inter- 
act harmonically with neighboring replicas, thereby form- 
ing a "chain" of "beads"—. The momentum distribution 
may be computed by means of an "open" chain, whereas 
position-dependent equilibrium properties are computed 
via "closed" chains. Although the vast majority of path 
integral simulations solely utilize "closed" paths, the mo- 
mentum distribution has been computed by means of 
open path integral simulation in superfluid heliumi^, and 
more recently in water-i2i2£iiSii22ii^. The latter results 
complement the measured experimental proton mome- 
tum distributions of water in the liquid and hexagonal 
solid phases*'^*'. 

The open path integral methodology has been recently 
implemented in conjunction with first principles molec- 
ular dynamics'^^ within the Car-Parrinello framework^'*. 
Previous computational studies of the momentum dis- 
tribution have not focused on the phenomena of proton 
tunneling or symmetric hydrogen bonds. First principles 
potentials facilitate the study of bond forming and break- 
ing events and thus make the study of these phenom- 
ena possible. Presently, we explore this topic by means 
of first principles open path integral molecular dynam- 
ics studies of high pressure ice. Although no neutron 
Compton scattering experiments have been performed 
on this system, there are several factors that make its 
study appealing. High pressure ice has been extensively 
studied via first principles molecular dynamics Simula- 
tions^'^e^'SS^^^'^^^^. In this work, we will con- 
sider three phases of high pressure ice. Ice VII, Ice VIII, 
and Ice X. By varying the volume of the simulation cells, 
interconversion of the phases may be induced. It has 



been previously shown in the pioneering work of Benoit 
and Marx'^'iSi that quantum delocalization and tunnel- 
ing play a crucial role in the transition between these 
phases. Each phase embodies a different hydrogen bond- 
ing "state," namely the typical case where the proton 
experiences relatively strong covalent and weak hydro- 
gen bonding (Ice VIII), proton tunneling along the hy- 
drogen bond from potential well to well (Ice VII), and 
the "symmetric" hydrogen bond (Ice X). Therefore, its 
position space distributions are already well character- 
ized in the literature. Furthermore, since bulk water is 
particularly amenable to previously developed open path 
integral methodologies;^!^ the problem is computation- 
ally tractable and several phases may be readily studied. 

We find that the uncertainty relation between the po- 
sition and momentum distributions of the proton^^ is 
evident in the behavior of each system. In particular, 
the more delocalized protons of Ice VII and Ice X pos- 
sess more localized momentum distributions in the hy- 
drogen bonding direction. In Ice X, the symmetric hy- 
drogen bond yields signatures of the momentum distri- 
bution that resembles those observed experimentally in 
other symmetrically hydrogen bonded systems^iii^. In 
the tunneling case, we observe a narrowed momentum 
distribution with an anomalous shape, namely one that 
is narrowed at low-momentum, but features tail behavior 
similar to the covalently bonded state. However, we find 
no clear indication of secondary peaks or features that 
have been associated with tunneling in a variety of ex- 
periments^ii^ii^iii^. Additionally, we present the resultant 
mometum distributions of a particle in a one-dimensional 
double well potential that show that secondary features 
may indeed be observed when tunneling is ground-state 
dominated. 

This article is organized as follows. In Section |lT] we 
review the methodology of first principles open path in- 
tegral molecular dynamics. In Section [IIII we discuss the 
three high pressure ice systems under study and the sim- 
ulation details, and then present the results of these com- 
putations in Section llVl An open path integral molecular 
dynamics simulation of a simple one-dimensional model 
of a particle in a double-well potential is presented in 
Section |V] and discussion and conclusions are given in 
Section IVll 

II. OPEN PATH INTEGRAL 
CAR-PARRINELLO MOLECULAR DYNAMICS 

The momentum distribution, n(p), may be written as 
the Fourier transform of the density matrix in position 
space, p(r,r'): 

n{p) = J drdr'e-^P ('-'''V(r,r') (1) 

The path integral discretization of the density matrix 
maps the quantum system onto a set of P replicas 
("beads") that obey classical physics, thereby allowing 
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one to utilize the machinery of computational classical 
statistical mechanics, namely Monte Carlo and molec- 
ular dynamics. The discretized density matrix may be 
written as: 



p{r,r') — lim 



dra . . . drp e-^^'^" 



(2) 



n—r 



2P 



p 



(3) 



This expression holds for a single particle but can be 
easily extended to multi-particle systems. 

Typically, one is interested in computing equilibrium 
averages of position-dependent properties. In this case, 
only diagonal elements of the density matrix are required 
and "closed" paths (for which r = r') must be sampled. 
However, the Fourier transform shown in Equation [1] re- 
quires off-diagonal components, and may be computed 
via "open" paths. In this formalism, the momentum dis- 
tribution is related to the Fourier transform of the end- 
to-end distance distribution of the open path. 

In order to improve the efficiency of sampling, closed 
path integral molecular dynamics simulations may em- 
ploy a coordinate transformation on {r} that decou- 
ples the harmonic interactions in the first term of Equa- 
tion [21 One such approach is the staging transforma- 
tion^i22i^i2S. The staging transformation has been re- 
cently extended to open path integral simulation^^. First, 
the following identity is employed: 



p ^ 

i=2 



X e 2^-^/3 



^ki-r-p+l 



where: 



[i - l)ri+i + ri 



rrii — m 



(4) 

(5) 
(6) 



and 



are known as the staging masses. We can 



then transform to the following set of coordinates, 
{ri,M2 . . .Mp+i, . . . ,rp+i} with Mi = r—r* fori = 2, P. 
The two endpoints of the chain undergo a transformation 
into relative and center-of-mass coordinates: 



Ml 

Up+l 



ri + rp+i 



2 

- rp^ 



(7) 

(8) 



The effective potential (Equation [S]) may be written as: 



2 >r;^ mfP 



2^12^2 



V{r,{{u})) + V{rp+,{{u})) 



E 



2P 

v{n{{u})) 
p 



where: {u\ 



i=2 

{m1,M2 . . .Mp,Mp+i} 



(9) 
(10) 



In order to generate molecular dynamics trajectories the 
forces fu — are computed as follows: 

(11) 
(12) 

1-2" 



fui 



fA + fB 
fA - fs 



pfr 



1 



ii = 2,P) (13) 



where: 



fE 



1 ^ \^ f^^^\ f 

-fsfri + "75 E 



i=2 



(14) 
(15) 



Staging open path integral molecular dynamics has 
been employed in conjunction with the Car-Parrinello^i 
scheme in Reference 1 2 2l This is a straightfowared exten- 
sion of pre-existent closed path integral Car-Parrinello 
molecular dynamics^SiiLii^. The corresponding extended 
Lagrangian is given by: 




V — 



gP 



dr 



X^mfP 2 

E^i«^i 



drcj,f'\r)4\r) - d,t 



with: 



9 = 




1 otherwise 



(16) 



(17) 



where /i is the fictitious electron mass, E is the ground 
state potential energy, and the final term enforces the 
condition of orthonormality upon the orbitals. The 
masses, m^, are associated with the velocities of the stag- 
ing coordinates and are chosen to be a multiple of their 
corresponding staging masses, mf'. 
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The momentum distribution is a single-particle prop- 
erty. For multi-particle systems within the open path 
integral formalism, it may be exactly computed when 
one path is opened and all others are closed. This re- 
quirement would lead to very inefficient sampling in bulk 
materials. However, it has been shown that if the paths 
of multiple particles are "opened" and these paths are 
sufficiently far apart from each other, the impact upon 
the resultant distribution is negligible^ ^. For a system 
of water molecules, one hydrogen per molecule is treated 
with an open path. All oxygen and all other hydrogen 
paths are closed. The end-to-end path distribution is 
then averaged over all open paths. 



III. SIMULATION DETAILS 

Ice possesses a rich phase diagram. At ambient pres- 
sure and below 0°C, water is most stable in a hexago- 
nal crystal structure. This is the form of ice for which 
the momentum distribution has been previously stud- 
^g j20, 21^22, 23 ^ However, under conditions of very high 
pressure, individual water molecules are arranged in in- 
terpenetrating cubic hydrogen bonded lattices. This ar- 
rangement forms an effective body centered cubic (BCC) 
lattice structure. In this study, we will concetrate on 
three phases. Ice VII, VIII, and X. Ice VIII is proton or- 
dered, exhibiting an anti-ferroelectric hydrogen bonding 
pattern. In comparison. Ice VII is proton disordered. Un- 
der higher pressures, the oxygen-oxygen distance reduces 
to the point where the proton's most stable position is 
equidistant between oxygen atoms and is located at the 
midpoint of the hydrogen bond axis. This "symmetric" 
form of ice is known as Ice 

The work of Benoit and Marx-S^j--^ has shown that by 
varying the lattice parameter (which changes the volume, 
and is equivalent to a change in pressure) of an Ice VIII 
cell one may, after a suitable equilibration period, gener- 
ate Ice VII and Ice X. In the case of Ice VII, the system 
will tunnel through the barrier along the hydrogen bond 
axis, thereby disrupting the proton ordering in the sys- 
tem. At even smaller volumes. Ice X becomes thermody- 
namically favored. A schematic is provided in Figure [1] 
that illustrates these concepts. 

Presently, we consider a 2 x 2 x 2 BCC supercell con- 
taining 16 water molecules at three different lattice con- 
stants. The lattice constants, as well as the correspond- 
ing molar volumes, pressures, and most probable oxygen- 
oxygen nearest neighbor distance are given in Table [H 
The approximate pressures are garnered from the equa- 
tion of state given by Hemley et al^. 

The first principles open path methodology is em- 
ployed in order to generate the trajectories. After an 
equilibration of 4 ps, each system is simulated for 75 
ps, with the exception of System 2, which is sampled 
for 120 ps. A timestep of 0.0725 fs is employed in all 
simulations. Each system is sampled at lOOK. The tem- 
perature is controlled by means of massive Nose-Hoover 



chain thermostat o ' i . Each path contains 32 repli- 
cas. The electronic states are evolved utilizing the Car- 
Parrinello methodology^! with a fictitous mass of 340 
atomic units. The electronic structure is described via 
the Kohn-Sham formulation of Density Functional The- 
ory^ where exchange and correlation effects are treated 
by the BLYP functional^i^. The valence orbital are 
expanded in a plane wave basis set with a cutoff of 75 
Rydberg. TrouUier-Martins norm-conserving pseudopo- 
tentials^ are utilized to model the valence effects of core 
electrons. The dynamical masses associated with the 
staging coordinates are set to be a factor of 4 larger than 
the staging masses. 

Despite the small number of water molecules in the 
simulation cell, there are 2048 electron states (32 repli- 
cas X 16 molecules x 4 states per molecule) present in the 
system. This is a relatively large system by the standards 
of first principles simulation, and only state-of-the-art 
computational resources make possible the calculation of 
the relatively long trajectories and multiple systems re- 
ported in this study. All computations are performed on 
IBM Blue Gene/L hardware with the CPMD program^—, 
which has been optimized for this archetecture^i^. 



System 
Number 


Lattice 
Constant 


Molar 
Volume 


Approx. 
Pressure 




1 


2.67 A 


5.74 cm^/mol 


90 GPa 


2.31A 


2 


2.84 A 


6.90 cm^/mol 


45 GPa 


2.45A 


3 


2.94 A 


7.62 cm^/mol 


31 GPa 


2.53A 



TABLE I: Characteristic values that relay the size of each 16 
molecule high pressure ice cell are given in the table above. 
The pressure is approximated from the equation of state given 
by Hemley et alM> The value of d^^ is the most probable 
oxygen-oxygen distance between nearest neighbor, hydrogen 
bonded molecules. 



IV. HIGH PRESSURE ICE 

The distributions in position and momentum space are 
computed in each system. As noted in Section |TIJ open 
paths are utilized for the computation of the momen- 
tum distribution, and closed paths are appropriate for 
the position distribution. Since our simulation contains 
both open and closed paths, we are able to use the closed 
paths for position distributions, and the open paths for 
the computation of the momentum distribution. The na- 
ture of this system in position space has already been 
elucidated in previous studies^'^'^. Here we repeat this 
work in order to explore the relation between the position 
and momentum space distributions. 

In Figure [H the first peak of the oxygen-oxygen ra- 
dial distribution function is shown. Shortening of the 
oxygen-oxygen distance is apparent as the molar volume 
is decreased. The position of the first peak of each distri- 
bution is reported in Table HI Although there is roughly 
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i 

icex 

FIG. 1: A schematic of the atoms involved in a single hy- 
drogen bond in the three high pressure ice phases presently 
under study. The gray circles represent oxygen atoms and 
the white circles represent hydrogen. As the pressure upon 
the system increases the average oxygen-oxygen distance de- 
creases, which has important consequences for the state of 
the proton. This may be covalently bonded (Ice VIII) , tunnel 
between wells (Ice VII) or lie in a symmetric state between 
the oxygen atoms (Ice X). 




distance (A) 

FIG. 2: (Color online) The first peak of the oxygen-oxygen 
radial distribution function in System 1 (solid curve) , System 
2 (dot-dashed curve) and System 3 (dashed curve). As one 
would expect, as the molar volume is decreased, the nearest 
neighbor oxygen-oxygen distance is as well. 



in a single well. This calculation showed that thermal 
hopping over the barrier is disfavored and quantum tun- 
neling dominates. As noted in Section [llll the tunneling 
disrupts the anti-ferroelectric ordering and engenders the 
formation of Ice VII. We note that the bimodal distribu- 
tion in Figure m is not perfectly symmetric. This may be 



two-tenths of an angstrom difference between oxygen- 
oxygen distances of Systems 1 and 3, this has a dramatic 
impact upon the nature of the proton that is confined 
on the potential energy surface. It is this shortening that 
drives the phase transition between the forms of ice under 
study.-^'^'^ 

The position space distribution of the proton along the 
oxygen-oxygen hydrogen bond axis is illustrated by the 
oxygen- hydrogen radial distribution functions (Figure [3]) 
and the probability distribution of the proton position 
along the hydrogen bond axis (Figured]). It can be seen 
that the proton in System 3 remains covalently bonded 
to its oxygen, although the covalent bond distribution 
is broader than in typical water phases. This system 
retains the Ice VIII structure. It can be seen in Figure 
[3] that the covalent bond and hydrogen bond peaks of 
the radial distribution function of System 1 merge. This 
broad single peak located at the midpoint between the 
two oxygen atoms is indicative of a symmetric hydrogen 
bond as found in the Ice X phase. 

Evidence of quantum tunneling can be seen in Sys- 
tem 2. The bimodal nature of the proton distribution in 
Figure m as well as the fact that one peak is near the co- 
valently bonded peak of System 3 indicates that there are 
tunneling events from one well to another along the hy- 
drogen bond axis. It was shown in the work of Benoit and 
Marx^°'^^ that classical protons at this molar volume and 
temperature do not cross the barrier and remain trapped 
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distance (A) 

FIG. 3; (Color online) The oxygen-hydrogen radial distribu- 
tion function in System 1 (solid curve), System 2 (dot-dashed 
curve) and System 3 (dashed curve). Whereas in System 3 
there is a distinction between covalent and hydrogen bonding 
distances, the two peaks have merged in System 1. 



We note that the present distributions are somewhat 
more delocalized when compared with the work of Benoit 
and Marxi^i^. This is likely a result of the use of a larger 
number of replicas in the present computation. However 
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FIG. 4; (Color online) The distance distribution of the proton 
along the oxygen-oxygen direction in System 1 (solid curve), 
System 2 (dot-dashed curve) and System 3 (dashed curve). 
This direction is analogous to the hydrogen bonding axis. 
One may note that the distribution of System 2 is delocal- 
ized across two wells. 



there are many other differences in the details of the sim- 
ulation that may impact this result, including trajectory 
length and the choice of exchange-correlation functional. 
Overall however, the description of the proton in position 
space along the hydrogen bond axis is in good agreement 
with this and later work.— 

The momentum distribution is plotted along the 
oxygen-oxygen axis (Figure O, as well along the two cor- 
responding perpendicular axes (Figure [6]) . These are ef- 
fective one-dimensional plots that are computed via the 
Fourier transform of the path end-to-end distance distri- 
bution along these directions. In Figure [SI one can view 
a trend that the momentum distributions in the direc- 
tions perpendicular to the hydrogen bond broaden with 
decreasing system molar volume. This is consistent with 
the uncertainty principle given that as the protons be- 
come more confined in position space, the corresponding 
momentum distributions have a greater variance. Aside 
from this difference, there is little distinction between the 
systems under study in these directions when compared 
to the momentum distribution projected onto the hydro- 
gen bonding axis (see Figure [5]). This is a logical conclu- 
sion as the large qualitative differences in position space 
occur in the hydrogen bond direction (see Figure [4|) , as 
shown presently and in previous work on high pressure 

j(,g30j31|35 

One may also note in Figure [5] that the distributions 
are similar along the two directions perpendicular to the 
hydrogen bond axis. This chemically intuitive result is 
in agreement with a previous study of the "shape" of the 
proton high pressure ice phases^, where it was found 
that the position space distribution in the perpendicular 
directions were of similar extent. In addition, the pro- 



ton's variance in the perpedicular directions was shown 
to decrease with increasing pressure"^^, thereby provid- 
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System 2 
System 3 
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momentum (inverse A) 



FIG. 5: (Color online) The proton momentum distribution in 
the oxygen-oxygen (00) direction in System 1 (solid curve), 
System 2 (dot-dashed curve) and System 3 (dashed curve). 
It is in this orientation that the distinctions between phases 
occur. 
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FIG. 6: (Color online) The proton momentum distribution 
perpendicular to the oxygen-oxygen direction (denoted "x") 
in System 1 (solid curve). System 2 (dot-dashed curve) and 
System 3 (dashed curve). Also plotted are the proton mo- 
mentum distributions in the mutually orthogonal direction 
(denoted "y") in System 1 (triangles pointing downward). 
System 2 (triangles pointing upward) and System 3 (circles). 
The differences in widths of these curves indicate the relative 
pressure upon each system. 



The position space distributions show that System 1 
contains symmetric hydrogen bonds, System 2 exhibits a 
bimodal proton distribution and in System 3, the protons 
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are covalently bonded. In Figure [5l we present the mo- 
mentum distributions in the hydrogen bonding direction. 
The covalently bonded System 3 possesses the narrow- 
est position distribution (see Figure 2]) and therefore the 
correspondingly broadest momentum distribution. The 
high-momentum tail of this distribution is dominated by 
the OH stretch (see Section |T|. In the symmetric hydro- 
gen bonded case (System 1), the more delocalized bond 
yields a narrower momentum distribution with a short- 
ened tail. This signature in proton momentum distribu- 
tions corresponds to a red-shift of the OH stretching fre- 
quency in stronger hydrogen bonded environments. This 
has been observed previously in the experimental* and 
simulation momentum distribution of liquid water and 
hexagonal ice^ . The symmetric hydrogen bond may 
be considered the "strongest" class of hydrogen bonding. 
Such an interpretation is bourne out by experiments on 
symmetric hydrogen bonds observed in water confined 
in nanotubes^^ and Rb3H(S04)2^® that exhibit greatly 
narrowed momentum distributions with shortened tails. 

The shape of the proton momentum distribution in the 
tunneling direction in System 2 lends to a more detailed 
description. It appears to have an anomalous shape when 
compared to the other distributions. Namely, it is narrow 
at low momentum, yet its tail behavior is similar to that 
of the covalently bonded System 3. This tail behavior 
is likely engendered by the localization in the covalently 
bonded well that is a component of the tunneling sys- 
tem. Therefore the highest frequency components of the 
system are similar to those exhibited in System 3. The 
narrowness exhibited in the low-momentum region is re- 
lated to the overall delocalized nature of the proton. The 
tunneling momentum distribution will be further inves- 
tigated in Section IVl 

In Figure [71 the spherically averaged momentum dis- 
tribution n[p), and the radial momentum distribution, 
p^n{p) are depicted for Systems 2 and 3. Note that the 
difference between the distributions is dominated by dis- 
tinctions present in the hydrogen-bonding direction, al- 
though this is somewhat muted by the contributions from 
the perpendicular orientations. 

The change in the spherically averaged momentum dis- 
tribution of Figure [7] reflects the transition in the state of 
the proton from System 2 (tunneling) to System 3 (co- 
valently bonded). Interestingly, this change bears sim- 
ilarity to that observed in the measured spherical mo- 
mentum distribution of a proton in the hydration shell 
of the globular protein lysozyme^^ at different temper- 
atures. Although these are different systems, the mo- 
mentum distribution is mostly dependent upon the local 
environment of the proton. Therefore it is likely that 
qualitative features are common among the set of pro- 
ton tunneling systems. To facilitate the comparison, we 
report the experimental distributions in Figure [51 The 
similarity between theory and experiment supports the 
interpretation that the observed change in the momen- 
tum distribution signals the onset of tunneling behavior. 
In our simulation, this is caused by a change of pressure 



whereas in the experiment, it results from a change in 
temperature. In both cases (see the insets of Figures [7l 
and [8]) the momentum distribution of a proton in a single 
well differs from that in a double well potential via a nar- 
rowed distribution in the low-momentum region. In ad- 
dition, the experimental tunneling distribution displays 
a distinct feature in the tail that is only visible when the 
distribution is multiplied by a factor of . This effect 
cannot be detected in our simulation. This may reflect 
genuine differences in the state of the proton in the ex- 
periment and in the present simulation. However, one 
should note that such small differences in the tail of the 
momentum distribution are beyond the precision of the 
current computation. 

In addition to the work of Senesi et al}'^ , experimental 
studies have reported secondary features in the momen- 
tum distribution-'^* or radial momentum distribution^>i^ 
in a variety of other systems, where they have also been 
interpreted as a signature of tunneling. The presence of 
secondary features in tunneling systems will be studied 
further in the next section. 

Experimental results such as those shown in Figure 
[51 manifest excess kinetic energy in comparison to liq- 
uid water and hexagonal ice^ i^^'^'^ . This is particularly 
expressed in broadened tails of the momentum distribu- 
tion. In some systems, this phenomenon has been identi- 
fied with tunneling modes®'^^. In contrast, excess kinetic 
energy is not apparent in the simulation results that are 
depicted in Figure [71 Instead, the tail of the distribu- 
tion in Systems 2 and 3 remains dominated by covalent 
bond stretching frequencies like those exhibited in the 
liquid and hexagonal crystal phases^ i^"'^^ . However, we 
note that from the simulation perspective, the kinetic en- 
ergy may be most accurately computed from closed path 
integral simulations^, as the tail of the momentum dis- 
tribution is difficult to obtain to high precision from the 
present methodology. 

V. A SIMPLE MODEL FOR TUNNELING 

In order to further investigate the tunneling behav- 
ior of the proton, we present an analysis of a particle 
in a one-dimensional potential. This degree of freedom 
corresponds to the displacement of the proton along the 
hydrogen bond axis as plotted in Figures [H and [5l Such 
models are often employed in the literature^, although 
we do not claim all the complexities of the high pressure 
ice system can be reduced to an effective one-dimensional 
form. Instead, it is simply a tool to easily investigate the 
equilibrium distributions of position and momentum in 
the tunneling regime. Therefore, we limit our discussion 
to an analysis of these properties and do not consider 
tunneling kinetics or dynamics. These issues have been 
extensively studied in the literaturci^iSi^'^i^i^i^ 

We utilize a potential of the following form: 

Viz) = ^mco^z^ + Ae-"'^'' (18) 
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FIG. 7: (Color online) The spherically averaged momentum 
distributions and radial momentum distributions of Systems 2 
(solid curve) and 3 (dashed curve) are plotted in panel (a) and 
panel (b), respectively. Each curve is normalized such that 
area under 47rp^n(p) is equal to one. The difference between 
the plotted momentum distributions is depicted in the inset 
of panel (a). 




with m = 1836, oj = 0.005, A = 0.012, and ^ = 0.0087. 
The parameter lo characterizes the confinement of parti- 
cle in the absence of the barrier, and A and ^ define the 
height and width of the potential barrier, respectively. 
All parameters are reported in atomic units and were 
chosen to yield a similar distribution to System 2 in po- 
sition space (see Figures HI and ITU)) . Our potential choice 
in Equation (TH] is by no means unique, as other model 
potentials may be employed^. 



FIG. 8: (Color online) The experimental proton momentum 
distribution of hydrated lysozyme. These results were re- 
ported in Figure 3 of Senesi et ali^ The spherically averaged 
and radial momentum distributions of a proton in a single 
well (dashed curve) and double well potential (solid curve) 
are plotted in the upper and lower panel, respectively. The 
normalization matches that of Figure[7l The curves above and 
below the distributions delineate the error in the measurement 
of n(p)^. The difference between the plotted momentum dis- 
tributions is depicted in the inset of the upper panel. Figure 
courtesy of R. Senesi. 



Given that the system is a single particle in a one- 
dimensional potential, we may directly diagonalize the 
system Hamiltonian in order to obtain the spectra. There 
are a large variety of schemes available that numerically 
solve the one-dimensional Schrodinger equation*'-. We 
have utilized SLEDGEj^ for the spectral computations pre- 
sented here. The shape of the potential barrier is plotted 
against the lowest eleven energy levels in Figure [HI It can 
be seen that the ground state and the first excited state 
lie below the barrier. Higher energy levels are nearly 
equally spaced as the harmonic term of Equation [TH] be- 
gins to dominate. 

The distributions in position and momentum space 
may be computed as a summation over the density ma- 



trix: 

n{x) = -Y^e-P^^Wdx)? (19) 

z^O 
^ oo 

n{p) = -Y.e-^^'\^M? (20) 

where the partition function is Z = X^i^o ■ This in- 

finite summation is truncated when e~^^^^^^^'' < 0.001. 
The populations of the lowest two energy states as well 
as the population of states above the barrier are given in 
Table [iTl The position and momentum distributions at 
various temperatures are depicted in Figure 1101 

Upon inspection of the position distributions, one 
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FIG. 9: (Color online) The model double well potential (solid 
curve) is plotted alongside the first eleven energy levels ol the 
spectra (dashed lines). Note that the ground and first excited 
state lie below the barrier. 



notes that the bimodal shape indicates there are barrier 
crossings at all temperatures. At low temperatures, only 
states below the barrier are populated, and quantum tun- 
neling is responsible for the traversing of the barrier. As 
the temperature increases, thermally activated crossings 
begin to contribute and there is a crossover from quantum 
to classical behavior^^. This crossover is not apparently 
detectable from the position distribution unless further 
simulations are performed on classical particles^° or iso- 
tope effects are studied^^. 

The momentum distribution provides complementary 
information to the position distribution and about the 
impact of quantum tunneling upon the system. In Fig- 
ure [TOl one finds that similar position space distributions 
may lead to rather different momentum distributions. A 
node is present in the momentum distribution at very low 
temperatures where the ground state dominates®^. These 
distributions are qualitatively similar to those obtained 
for KDP in neutron Compton scattering experiments^^, 
as well as those garnered from ground state models that 
are described in the literatur e^'^^i^'^ . The node disap- 
pears at intermediate temperatures where a non-gaussian 
momentum distribution with an extended tail is present. 
This is the general shape of the momentum distribution 
of System 2 (see Figure [S]) . At higher temperatures the 
distribution becomes more gaussian-shaped in the regime 
of the quantum/classical crossover. 

In addition, we note that the presence of nodes in- 
dicates a barrier in the underlying potential energy sur- 
face but not necessarily quantum tunneling of the ground 
state. A perturbation analysis predicts a node in the 
ground state for any barrier of finite height, even in the 
case where the barrier is "washed out" by zero-point- 
motion and a unimodal position distribution is present 
(see Appendix [A| . However, for systems that exhibit bi- 
modal position distributions, a node in the momentum 



distribution is a signature of ground ("coherent") state 
tunneling. Mixed ("incoherent") state tunneling is indi- 
cated by a momentum distribution of shape similar to 
that in the model at T = 300K (Figure [TH]) or in System 
2 (Figure p. 



r 


30K 


lOOK 


BOOK 


lOOOK 


2400K 


Po 

Pi 

P2+ 


100% 
0.00% 
0.00% 


99.1% 
0.881% 
0.00% 


82.8% 
17.2% 
0.00% 


57.1% 
35.6% 
7.29% 


37.5% 
30.8% 
31.7% 



TABLE 11: This table depicts the populations of the ground 
state, first excited state, and all other levels as a function 
of temperature for a single particle in the potential given in 
Equation [181 
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FIG. 10: (Color online) The position (left panel) and mo- 
mentum (right panel) distributions of the double well model 
at several temperatures. One may note that there are larger 
qualitative difi'erences in the momentum space picture than 
the position space picture with the increase in temperature. 
Note the disappearance of secondary features in the momen- 
tum distribution which is evident at T = 300K. 

To further investigate the temperature dependence of 
the model system and its relation to the high pressure 
ice calculation, open and closed path integral molecu- 
lar dynamics simulations were undertaken at T=30K, 
T=100K, and T=300K. Note that at each of these tem- 
peratures, there is virtually no population of energy lev- 
els above the barrier and tunneling effects dominate. 
The position and momentum distributions garnered from 
these computations are in good agreement with those 
presented in Figure 1101 The path may be characterized 
by the root mean square imaginary time correlation func- 

7^2(T) = ^[z(t)-z(o)]^^ o<T<f]n (21) 

The time independence ("flatness") of the root mean 
square imaginary time correlation function near the mid- 



point is related to the degree of ground state domi- 
nance^^. In Figure [TT] we plot this function for the three 
simulated temperatures and System 2. The T=30K curve 
is the flattest, indicating that the system is in the ground 
state, whereas the others show some time dependence 
that is suggestive of the presence of some excited state 
wei 
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FIG. 12: (Color online) The distribution of the centroid of 
the paths (left panel) and the distribution of the radius of 
gyration (right panel) of the double well model at T=30K 
(solid curve), T=100K (dashed curve) and T=300K (dotted 
curve with circles) plotted against that of System 2 in the 
hydrogen bonding direction (dot-dashed curve). 



FIG. 11: (Color online) The root mean square imaginary 
time correlation function of the model system at T=30K (dot- 
ted curve with circles), T=100K (dotted curve with squares), 
T=300K (dotted curve with diamonds), and System 2 of the 
high pressure ice calculation (dotted curve with triangles). 
The imaginary time is reported in units of /3fi. 



In addition to indicating ground state dominance, the 
height of root mean square imaginary time correlation 
functions is related to the degree of path localization. In 
order to further characterize this property, the distribu- 
tion of the centroid and the radius of gyratio n^^'''" of the 
proton paths is plotted in Figure [T^l The centroid is 
defined as the center of mass of the path: 

1 P 

and the radius of gyration is defined as: 

1 ^ 

1=1 

In this way the average and variance of the path may be 
characterized. 

In the plots of the radius of gyration and the cen- 
troid, we see qualitatively distinct behavior in the one- 
dimensional model at each temperature. When the sys- 
tem is in the ground state (T=30K), the broad paths are 
centered about the midpoint of the position distribution 
and delocalized across the two wells. This finding is in 
agreement with Feynman-Kleinert perturbation analysis 
of a double well potential that display a minimum in the 
centroid potential of mean force at the barrier at zero 
temperature^!. In the high temperature case, the spread 



of the path has greatly narrowed, and its centroid distri- 
bution is now bimodal, centered about the two wells of 
attraction. Therefore, the tunneling appears to be occur- 
ring by two different "mechanisms" in the ground (30K) 
and mixed states (BOOK), a single delocalized species in 
the former and "path hopping" in the latter. At lOOK, 
the system is in an intermediate state between these two 
situations, as can be starkly seen in the bimodal form 
of its radius of gyration distribution (see Figure [T2| . 
The path hopping at higher temperatures is made fa- 
vorable by the stiffened harmonic interactions between 
beads at higher temperatures that force the localization 
of the path. In the high temperature (classical) limit, one 
would expect the path to collapse to a point and ther- 
mally hop between wells. The behavior at 300K seems to 
be approaching this limit, although as remarked earlier, 
quantum tunneling still dominates. Delocalized paths 
centered at the top of the barrier and localized paths in 
the potential wells of a double well potential were char- 
acterized in the work of Tuckerman and Marx^*, and this 
is in good agreement with the picture presented here. 

The centroid and radius of gyration distributions of 
System 2 along the hydrogen bonded axis is plotted in 
Figure[T21 One may see that although this simulation was 
performed at lOOK, the results are qualitatively similar 
to those at 300K. Both System 2 and the BOOK model mo- 
mentum distributions lack secondary features and have 
similar shapes (see Figures [51 and [TU)) . The similarity in 
position and momentum space with the model at BOOK 
suggests that System 2 is in a mixed state tunneling 
regime. We remark that although the model potential 
of Equation 1181 was tuned to mimic System 2, there are 
still effects that cannot be captured by the model, such 
as the proton correlations that arise from the enforce- 
ment of ice rule. Despite these limitations, the model 
potential provides a very useful qualitative picture of the 
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temperature dependence of tunneling phenomena. The 
"crossover" from ground state to mixed state tunnehng, 
as indicated by the population of the energy levels, is 
apparent in the shape of the momentum distribution, as 
well as the spread and centroid of the particle's path in- 
tegral representation. 

Finally, we remark that our one-dimensional model is 
characterized by a potential that does not depend on tem- 
perature. Within this model, tunneling features in the 
momentum distribution become more pronouced at lower 
temperatures. By contrast, experiments on systems such 
as KDP— and the hydration shell of lysozymei^ have 
reported secondary tunneling features that only exhibit 
themselves above a critical temperature. This behavior 
must be facilitated by a structural transition above the 
critical temperature in which the potential experienced 
by the proton, and particularly the tunneling barrier, 
changes significantly. Furthermore, we note that the ex- 
perimental interpretations rest in part on a ground state 
modeli^ii^, in agreement with the present finding that 
nodes in the momentum distribution are associated with 
ground state tunneling. 



VI. CONCLUSIONS 

In this work we have presented an investigation of the 
position and momentum space distributions of the proton 
in tunneling and symmetric hydrogen bonded systems. 
Novel first principles open path integral molecular dy- 
namics algorithms were utilized in order to compute the 
momentum distributions. Three phases of high pressure 
ice were studied at lOOK. Each phase typifies a quali- 
tatively different state of the proton, covalently bonded 
(Ice VIII), tunneling (Ice VII), and equally shared be- 
tween nearest-neighbor oxygens (Ice X). 

In Ice X, symmetric hydrogen bonds are present. This 
phase is characterized by a broadened position space 
when compared to the covalently bonded Ice VIII dis- 
tribution. In accordance with the uncertainty relation 
between position and momentum, the computed momen- 
tum distribution is narrowed in the direction of the hy- 
drogen bond axis. This is also in agreement with the 
experimental signatures of symmetric hydrogen bonding 
that have been observed in other systemsiiii^. Further- 
more, this phenomenon is an extreme case of the narrow- 
ing of the proton momentum distribution that accompa- 
nies the red-shift in the oxygen-hydrogen stretching fre- 
quency of hexagonal ice relative to liquid wateri^i2^. 

The resultant proton momentum distribution of the 
tunneling system (Ice VII) also shows some qualitative 
correspondence to what has been observed experimen- 
tally. The tunneling distribution is narrowed when com- 
pared to the covalently bonded distribution, in accor- 
dance with the relatively delocalized bimodal distribution 
that is a signature of tunneling in position space. The 
tail behavior is similar to that of the covalently bonded 
distribution, indicating that the high frequency associ- 



ated with the covalent bond is still present in the system. 
However, we did not observe tunneling signatures such as 
nodes or secondary features in the momentum distribu- 
tion that have been dectected in experiment. As noted 
in Section llVi such features may be rather subtlc^"^ and 
require sensitivity in the tail of the momentum distri- 
bution where the intensity is very small. Such precision 
is beyond the scope of the currently utilized methodol- 
ogy, and addressing this issue is an important goal of 
future algorithmic development. Furthermore, such fea- 
tures intimately depend upon the details of the potential 
experienced by the proton. Direct comparison of theory 
and experiment for a system in which these signatures 
are unambiguously detected should therefore provide an 
extremely accurate test of the current description of hy- 
drogen bonding. 

To better understand how tunneling affects the mo- 
mentum distribution, we performed an open path inte- 
gral simulation of a particle confined in a one-dimensional 
double well potential. The position and momentum dis- 
tributions may be directly computed in this system from 
the spectrum of eigenvalues and eigenvectors. We found 
that a node is present in the momentum distribution at 
low temperature, which is indicative of ground ("coher- 
ent") state tunneling. In this case, the form of the dis- 
tribution is in qualitative agreement with experimental 
studies of potassium dihydrogen phosphate-^. The Ice 
VII momentum distribution bears a resemblance to the 
results of the one-dimensional model at higher temper- 
aure (300K), where the node is "washed out." A clue as 
to why is revealed upon inspection of the centroid and 
radius of gyration distributions of the path, which show 
that the path tends to be localized in the wells rather 
than delocalized across the domain, as path integral sim- 
ulations of the low-temperature (node-containing) one- 
dimensional systems. The relative localization in posi- 
tion space lends itself to a broader distribution that is 
somewhat closer to that of the covalently bonded sys- 
tem. 

Whereas in a one-dimensional double well potential a 
secondary feature is always present when coherent tun- 
neling dominates, in a three dimensional system such fea- 
tures may be washed out in the spherically averaged dis- 
tribution. This suggests that these characteristics may 
be more easily detected in crystalline systems where the 
momentum distribution can be measured along high sym- 
metry directions. On the other hand, the momentum 
distribution is very sensitive to the spatial environment 
experienced by the protons and simple one-dimensional 
models, while useful for identifying tunneling features, 
may be too crude for realistic predictions. Indeed, the 
qualitative differences in the momentum distribution of 
the model and Ice VII at lOOK, while certainly reflecting 
a crude parameterization of the model, may also indicate 
that important features of the potential energy surface of 
the high pressure ice system are not reducible to a single 
particle, one-dimensional form. 

Finally, as noted in Section U and discussed in Sec- 
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tion IIVI some experiments report excess kinetic energy 
as compared to ambient liquid water~i^iH. In this work, 
such behavior is not apparent. AU calculated tunneling 
distributions possess tail behavior that is similar to the 
covalently bonded Ice VIII (System 3). The highest fre- 
quency in the system remains associated with the oxygen- 
hydrogen stretch, and no higher frequency modes act in 
the tunneling direction. This result should be contrasted 
with the interpretation of recent experimental data sug- 
gesting that excess kinetic energy is associated with tun- 
neling modes^iSiii^. Further investigations are necessary 
in order to unravel these findings. 



Acknowledgments 

We would like to acknowledge C. Andreani and R. Sen- 
esi for useful discussions and to thank R.S. for providing 
us with Figure [H J. A.M. acknowledges the Fannie and 
John Hertz Foundation for its support of his graduate 
work. Partial support for this work was provided by the 
DOE under grant DE-FG02-05ER46201 and by the NSF- 
MRSEC program through the Princeton Center for Com- 
plex Materials (PCCM), Grant DMR 0213706. In addi- 
tion, we would like to acknowledge Princeton University 
and IBM for the use of their computational resources. 

This article has been submitted to the Journal of 
Chemical Physics. After it is published, it will be found 
at _http : / /j cp . aip .org 



APPENDIX A: PERTURBATION ANALYSIS OF 
THE GROUND STATE MOMENTUM 
DISTRIBUTION 

In this Appendix, we utilize a perturbation analysis 
in order to show that the ground state momentum dis- 
tribution always possesses extra nodes when a potential 
barrier exists. This is equivalent to showing that the 
Fourier transform of the ground state has extra nodes. 

For example, this may be seen in Figure [13] where the 
momentum distribution of the model potential (Equation 
[T5)l with A — 0.004 and all other parameters unchanged 
from Section|V]is plotted. For these parameters, the bar- 
rier height is 300K and the ground state energy is 585K. 
Therefore, zero-point motion washes out the barrier and 
the position distribution is unimodal. However, a node 
is still detectable in the momentum distribution, though 
the secondary feature is diminished with respect to larger 
values of A. 

Starting from a harmonic oscillator 

V = imu;2^^ (Al) 

the eigenvalues are i?„ = {n + ^)uj, and corresponding 
eigenstates are denoted by \n). The ground state wave- 
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FIG. 13: (Color online) Ground state position (left panel) and 
momentum (right panel) distribution for the low barrier dou- 
ble well potential. The inset depicts that the secondary fea- 
ture persists in the momentum distribution even when zero- 
point motion washes out the barrier in position space. 

function in position space is: 

(z|0) = (^—'^ exp (-mcjzV2) Ho{^/^z), (A2) 

The functions iJ„ are Hermite polynomials. The Fourier 
transform of the ground state is: 

m-(-y^'J^e-^. (A3) 
V TT / V muj 

We are interested in the Fourier transform of the ground 
state in the presence of a potential perturbation: 

AV^Ae-""^''. (A4) 

For simplicity, m,Lo,^ are chosen as before, and only A 
is varied. 

If 1^1 is small, the change in wavefunction can be cap- 
tured by the first order perturbation. Denote by \lp) the 
new ground state, we have 

i even and i>0 

Note that since the external potential is an even function 
all the odd terms in the sum vanish. We note that the in- 
finite summation in Equation (jASP cannot be reduced to 
a simple form and is therefore difficult to use for further 
analysis. Despite the slow decay of this series, numer- 
ical investigation shows that qualitative information is 
already captured by including only the first two terms 
(|2) and |4)) in the first order expansion. In particular, 
the existence of the node in the presence of an arbitrary 
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small potential barrier is already discernible within this 
approximation. 
One finds: 



, , . , , ,1 /TOW \ 1/2 



(A6) 



V2 aV2 



Here a = 1 + — . 
Similarly: 



(A7) 



Also we have: 



, fmijj\ 



1/4 1 



V2(tow)3/- 



-e-^{muj -2p^) (A8) 



^4 = 

V TT 



1/4 



/ %/24(tow)S/2 



e 2" 



(A9) 



Use £"„ = (« + l/2)tj, and substitute all the above 
relations into the first order perturbation with the first 
two even terms, we have: 



1 



1/4 



e 2r, 



1 



A{a - 1)2 / p 



Here pi and p2 are defined as: 



Pi 



P 

niLO 



Pi 



P2 = 



V6q;2 + 4q + 6 - (q + 3) 
2(a- 1) ' 
-V6a2 +4a + 6 _ + 3) 
2(a- 1) ■ 



-)}■ 

(AlO) 



(All) 
(A12) 



Since a > 1 by definition, the relations above are always 
well defined. Also note that: 



Ga^ + 4a + 6 - (a + 3)^ = (a - l)(5a + 3) > 0, (A13) 



and it holds that: 



Pi > 0, P2< 0. 



(A14) 



The potential barrier corresponds to the case when 
A > 0. No matter how small A is, there always ex- 
ists a p^ /{rauj) > pi and Equation (IA10|) equals to zero. 
The position of the node can be easily calculated from 
Equation (jAlOp . For instance, for A = 0.004 (the low 
barrier case), Equation (jA10|) gives the node position at 
p = ±12.99. Compared to Fig. [HI Equation (^10)) gives 
a good prediction of the position of the extra node. 

It is also of interest to look at the case of A < 0, that is, 
when the potential barrier is substituted by a potential 
well. The factor inside the braces in Equation (|A10[) 
reaches its minimum at p = 0, and this minimum is 



A{a-1)^ 
1 :-^77^ P1P2 



(A15) 



Note that piP2 is a fixed number and only depends on 
a, therefore when A is small, this minimum is always 
positive. Then it follows that {p\ip) > everywhere, and 
there is no node. For example, when A = —0.004, the 
minimum is 0.94 and far from zero. 
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